Specular Electron Focusing between Gate-Defined Quantum Point Contacts in Bilayer Graphene

We report multiterminal measurements in a ballistic bilayer graphene (BLG) channel, where multiple spin- and valley-degenerate quantum point contacts (QPCs) are defined by electrostatic gating. By patterning QPCs of different shapes along different crystallographic directions, we study the effect of size quantization and trigonal warping on transverse electron focusing (TEF). Our TEF spectra show eight clear peaks with comparable amplitudes and weak signatures of quantum interference at the lowest temperature, indicating that reflections at the gate-defined edges are specular, and transport is phase coherent. The temperature dependence of the focusing signal shows that, despite the small gate-induced bandgaps in our sample (≲45 meV), several peaks are visible up to 100 K. The achievement of specular reflection, which is expected to preserve the pseudospin information of the electron jets, is promising for the realization of ballistic interconnects for new valleytronic devices.

bulk crystals [3] and picked up with a PC layer at temperatures between 60 and 90 • C. The resulting heterostructure was released on a clean SiO 2 substrate with Au markers by melting the PC layer at temperatures above 150 • C. Subsequently, the PC covering the stack was removed from the surface by dissolving it in chloroform. The stack was then annealed for 1 h in an Ar atmosphere at 400 • C before contact preparation. At this stage, an atomic force microscopy (AFM) image was taken in the AC mode of a Cypher AFM to determine the thickness of the different flakes (Fig. S1a) The contacts to the bilayer graphene were defined using e-beam lithography. After defining the contact pattern, the top hBN was etched with a mixture of CHF 3 and O 2 with a 10 to 1 flow ratio, 40 W of power and a pressure of 5 µbar [4]. This recipe gives an etch rate for hBN of approximately 30 nm/min and leaves the BLG edges exposed. After etching, the Ti(5 nm)/Au(35 nm) electrodes were deposited using e-beam evaporation. The top gates were prepared using the same method replacing the CHF 3   In each point the current is applied following the sequence: I 1 = +I, I 2 = −I, I 3 = +I while the nonlocal voltage is measured after waiting 100 ms at each stage. The nonlocal voltages used here are obtained using V anti = (( This technique allows us to correct for any background signals. For completeness, is also extracted from the measured data and can be retrieved from [5]. The TEF measurements at different T were performed using a standard lock-in technique at a frequency f ≈ 18 Hz.

OSCILLATIONS
To determine the capacitance of the backgate (C bg ) in an accurate way we have used Shubnikov-de Haas oscillations. In a four-terminal measurement configuration, we measured where e is the electron charge, h is the Plank constant and B i and B i+1 correspond to the adjacent field positions where R 4p is minimal, indicated as crosses in Fig. S2a. The result from Equation S1 is shown in Fig. S2c and plotted vs. V bg as blue dots. Because the carrier density changes with the gate voltage following n = C bg (V bg − V cnp )/e, where V cnp is the position of the charge neutrality point, we have used a linear fit n = AV bg + B (orange line) to extract C bg = Ae and V cnp = −B/A. With ϵ hBN = C bg t hBN /ϵ 0 , where ϵ 0 is the vacuum permittivity, and the bottom hBN thickness from AFM measurements (t hBN ≈ 23 nm), we estimate its dielectric constant, ϵ hBN ≈ 3.75, in agreement with [6].

S3. ELECTRONIC MOBILITY
To estimate the device mobility, we have measured the BLG channel resistance as a function of V bg . The measurement geometry is shown in Fig. S3a and the result for the square resistance R sq = V × W/(I × L) is shown in Fig. S3b. Here, W = 7.2 µm is the BLG width between the V probes and L = 2 µm is the separation between them. In addition to the expected trend with a peak at the charge neutrality point (CNP) [3], we observe a dip at V bg = −1 V. From the two-point measurements used to characterize the QPCs (which show a peak for V bg ≈ −1 V), we conclude that it corresponds to the CNP of the BLG near the contacts. Thus, we attribute the drop in the measured resistance to the gate-tunable invasiveness of the voltage probes, which is minimal for V bg ≈ −1 V, when the contact resistance is maximal, resulting in an enhancement of the effective mobility. This explanation is consistent with the channel being ballistic and the contacts having a small overlap with the current path.
We have also plotted σ = R sq −1 vs. V bg , which does not show the linear trend with V bg as expected from the Drude model for constant mobility (µ), that predicts σ = neµ. This shows that the effective µ depends on n. To determine µ (Fig. S3c) we have used µ = σ/(ne) and used µ to estimate the momentum scattering time τ p = m * µ/e, where m * = 0.034 × m e is the effective mass in BLG [7] and m e the electron mass. With v f ≈ 0.5 × 10 6 m/s, the meanfree-path l mfp = v f ×τ p ∼ 1 µm, is comparable to the 2 µm separation between the V probes.
This observation indicates that the channel is in the ballistic regime and the measured values represent a lower bound to the actual device quality. The underestimation of τ p using this method is confirmed by the clear observation of focusing between QPCs placed at a distance of L =4 µm, requiring a ballistic path of L × π ≈ 12.6 µm. Finally, there is not a significant difference between the electron and hole mobilities for |V bg | > 1 V, indicating that the electron-hole asymmetry in the focusing signals is not caused by a difference in mobility. Firstly, the collector resistance (R c ) decreases with increasing temperature. This decrease is caused by the small bandgap opening at the BLG under the split gates (around 45 meV for an applied electric field of approximately 0.5 V/nm [8]), that exhibits thermally activated behavior and has to be taken into account to analyze the T -dependence of the focusing signal. Secondly, the magnetoresistance below 0.5 T is smaller than 10% in both cases, allowing us to analyze the low-B focusing peaks assuming that R c is constant through the

B-sweep.
For V bg = −3 V, the low-T data shows a small peak at B = 0, that is absent for V bg = +3 V.
This peak resembles weak localization (WL). Its large width indicates that, if its origin is WL, it must come from a region with a small phase coherence length (of the 100 nm range).
We conclude that it most likely originates from WL near the contacts between the BLG and Ti/Au electrodes, which have been doped by the etching process and are expected to have worse transport properties than the rest of the channel.   Figure S4. Quantum point contact magnetoresistance at V bg −3 V (a) and V bg = +3 V (b).

A. Different temperatures
To determine how the scattering rate changes with T we have extracted the area under the focusing peaks shown in Fig. 4 of the main manuscript. This operation requires the assumption that the reference electrode, which is placed far from the injector, only contributes to a smooth background in the measured signal. Such assumption is justified considering the long distance between this electrode and the injector and the large effective width it has due to the absence of a pair of confining top gates.
The area under the peaks has been extracted following several steps: R nl has been normalized by R c .
A linear background has been corrected from the data.
The data has been slightly smoothed (see the small difference between the black curves and colored scatters in Fig. S5a-d).
The minima in the nonlocal signal have been identified using the find peaks function from the python package scipy.signal on the reversed data (-R nl /R c ).
The background has been defined for each peak by linear interpolation between the extreme points.
The area between the data and background has been calculated using the trapz function from the python package numpy.
The scattering rate τ −1 p has been calculated using [9]: where v F ≈ 9.0 × 10 6 m/s is the Fermi velocity of BLG at V bg = +3 V, A p (T ) is the area under the p-focusing peak (starting from B = 0) at each T and A p (T base ) is the area of this peak at T = 2 K.
The result of the process described above is shown in Fig A quadratic T -dependence of τ −1 p is associated with electron-electron interactions [9].
Thus, our analysis indicates that both scattering terms are relevant. By calculating T 0 = b/a, which is the T where the quadratic term starts to dominate over the linear term, we obtain  T is fit to aT 2 + bT + c, whereas in (g) and (h) it is fit to dT 2 . The fitting parameters can be found in Table S1. The inset of panel b shows the ratio Since we are analyzing the p = 2 peaks, it would be tempting to attribute the electronhole asymmetry to diffuse scattering at the edge (DSE), that the TEF spectra indicates may be stronger for holes. If DSE was T -dependent, it could lead to a faster p = 2 peak decay with T and artificially enhance t −1 p . By monitoring the T -dependence of A 2 /A 1 we can determine whether the bT term is dominated by T -dependent DSE because, in this We attribute it to the thermally-activated transport across the weakly-gaped (∼ 40 meV) BLG region at 100 K where k B T ≈ 8.6 meV.

Effect of the R c -dependence with temperature
The analysis shown above, which assumes that the τ −1 p = aT 2 term is caused by electronelectron interactions, is based on the results from Ref. [9], where the following result is obtained for monolayer graphene and a hard wall potential: where α ≈ 0.518, β ≈ 2.28, v = 1 × 10 6 m/s is the Fermi velocity of graphene, k B is the Boltzman constant, k F is the Fermi wavevector, and A = w 2 c + w 2 i , where w (i)c is the (injector) collector QPC width that determines the electron collimation. This implies that a change of the widths with increasing T could affect our analysis. To quantify this effect, we have assumed that the changes in R c with T shown in Fig. S4 are due to the activation of the BLG regions near the QPC which have a reduced bandgap. Under this assumption, the changes in R c should correspond to changes in w = w i ≈ w c following R c ∝ 1/w. If we assume that w = 50 nm, the separation between the split gates, and take an upper bound to the change in QPC resistance with T of 50%, we obtain a correction to the log in Equation S3 of 18% at 100 K. Because w is used to determine the electron jet collimation [9] and electron jets emitted by QPC have been shown to be collimated [12], this correction is expected to be an upper bound. However, to make our analysis more robust against this effect, we have performed fits to τ −1 p = aT 2 + bT + c for V bg = +3 V limiting the T range to 80 and 60 K, where the correction is expected to be of less than 5%. In the former case, T 0 = 50 ± 10 K has been obtained and, in the latter, T 0 = 70 ± 10 K. The fact that both values are below 100 K indicates that the τ −1 ee term dominates over the linear term at 100 K, consistently with Ref. [13].

B. Base temperature
To determine the specularity of the edge reflection in the TEF measurements, we have determined the area under the different peaks in the R nl data shown in Fig. 2c of the main manuscript. To take into account the contact magnetoresistance shown in Fig. S4 we have normalized R nl by R c using the data from Fig. S4 at T = 10 K. The result, obtained following the procedure described in Section S5 A, is shown in Fig. S6 and shows that, for V bg = −3 V, the peak amplitude decays much faster than for V bg = +3 V, as shown in the main manuscript. Additionally, the amplitude of the V bg = +3 V signal is still around 25% larger than the V bg = −3 V case. We attribute this small difference to the fact that the detector has a slightly larger asymmetry at V bg = ±3 V than the injector, which is the contact we corrected for. At the inset of Fig. S6b we show the normalized area under the different peaks in both V bg = −3 V and V bg = +3 V cases where one can see more clearly the faster decrease of peak amplitude in the former.
An additional feature which can be identified in the V bg = +3 V data is the apparent beating pattern which results in the splitting of the p = 4 TEF peak. To infer whether it is compatible with the expected interference pattern arising from TEF between QPCs which The result, which is compatible with the first six peaks, indicates that a beating pattern with two separate frequencies can explain most of the spectra obtained at V bg = +3 V.

S6. TRANSVERSE ELECTRON FOCUSING AT DIFFERENT GEOMETRIES
In the presence of trigonal warping, the transverse electron focusing (TEF) spectra are expected to depend on the relative orientation of the QPCs with respect to the crystallographic directions of the BLG [14,15]. For this reason, we have studied TEF using gatedefined QPCs which are oriented in different directions. In particular, the geometry used to obtain the TEF spectra shown in the main manuscript, which we also show in Fig. S7d, has a rotation of 30 • with respect to the second set of QPCs, which is shown in Fig. S7e (see Fig. S1c for an overview of the whole sample). Since the QPCs in area B (Fig. S1c) are oriented along a straight edge in the BLG, which is most likely a crystallographic direction [16], the QPCs in area A are expected to be aligned along an armchair (or zig-zag) direction while the ones in area B must be along a zig-zag (armchair) direction, although we cannot tell which one is which. To find if the focusing peaks occur at the expected B, we use [17]: where ℏ is the reduced Plank constant, k f the Fermi wavevector, θ the electron incidence angle from the QPC with respect to the normal, and L the contact separation. Note that Equation S4 corresponds to Equation 1 of the main manuscript.
To compare Equation S4 with the measured data, we have identified B f as the B values C2 C2L Figure S7. TEF measurements obtained using the QPCs in rectangles A and B of Fig. S1c. (a-c lines are the fits to B = B 0 + (dB/dp) × p, where B 0 accounts for the magnet remanence and dB/dp is the slope. (j-l) Slope of fits obtained from panels g-i, respectively, together with a fit to Equation S4 assuming normal incidence (θ = 0; dashed lines). where R nl is maximal, plotted B f vs. p at each V bg , and fitted the data to B f = B 0 +(dB/dp)× p, where B 0 accounts for the coercivity of the magnet and dB/dp quantifies the change in B f with p (see Figs. S7g-i). After this, we have plotted |dB/dp| vs. n and compared this result with Equation S4 assuming cos(θ) = 1, which provides the maximal peak separation (see Figs. S7j-l). We observe that the separation between the measured peaks is even larger than predicted by Equation S4. This result is more pronounced in Fig. S7k, and may be consistent with an impurity obstructing the charge transport path and leaving a smaller peak at B < B f (which we did not consider in Fig. S7h) and a larger peak at B > B f . As a result the spacing between the first and second peaks is artificially enhanced and, since the p = 3 peak is not clear enough, |dB/dp| is overestimated. The opposite occurs when we consider the smaller peaks at lower B (Fig. S7k, orange dots).
For V bg > 0 the background signals in Fig. S7e are much larger than in Fig. S7d and additional features that are not expected from TEF are observed. We attribute them to quantum interference and restrict our comparison to V bg = 3 V. In this case, up to eight TEF peaks are visible, as in Fig. S7d, but without the splitting of the p = 4 peak. Looking at the spacing between the peaks, we see that in both cases it is close to the fit for cos(θ) = 1, indicating that there is only one current jet which departs at normal incidence from the QPCs.
For completeness, we have also measured TEF over a distance of 4 µm using the same injector as in Fig. S7e but connecting the detector to the lowest QPC electrode. The result is shown in Fig. S7f and, as expected from Equation S4, the peak spacing is approximately halved with respect to Fig. S7e. The result is summarized by comparing Fig. S7k with Fig. S7l. We note, however, that in both cases the spacing obtained from the spectra is slightly larger than predicted by the model.
To study the influence of the size quantization of the QPCs on the TEF data we have also prepared some QPCs with a horn-like shape which do not show size quantization for any and S8d, inset, the electrodes used for TEF in the C3 and C4 geometris are also rotated by 30 • . We thus believe that the close similarity between the TEF spectra obtained in both cases indicates that trigonal warping does not play a dominant role in our measurements.

A. Reciprocity
To confirm that the TEF measurements are in the linear response regime we have measured R nl in C1 (R C1 nl ) and its reciprocal geometry (C1R, R C1R nl ), obtained by swapping the current and voltage leads [18] Figure S9. (a) and (b) Measurement geometries C1 and C1R used to measure the reciprocity of the TEF data, which is shown in panels (c) and (d) for V bg = +3 V and −3 V, respectively. The B has been reversed for the C1R measurements to show the agreement between both curves in a more clear way.

ITY POINT
When measuring R nl at |V bg | < 1 V the double-gated regions have a very small band gap and play a role in charge transport. We have measured R nl vs. B in configuration C1 of the main manuscript and the results are shown in Fig. S10. In Figs. S10a and S10b, V tg has been set to keep the double-gated regions at n ≈ 0. In Fig. S10c  map in Fig. 1b of the main manuscript). A similar behavior is observed in Fig. S10a. In this case, an additional effect occurs: R nl increases dramatically when |B| > 1.5 T. We believe that this is due to the formation of an insulating state in BLG near the CNP at high B [19].

S8. TRANSVERSE ELECTRON FOCUSING SAMPLE 2
We have also measured TEF in a second heterostructure. Its top and bottom hBN thicknesses are 49 and 40 nm, respectively, and its optical microscope image is shown in  applied current, defined in Fig. S11a) is shown in Fig. S11b as a function of V bg , applied to a multilayer graphene backgate and V tg , applied to the top gates surrounding the left QPC.
Following the protocol described in the main manuscript, we extracted G of the left and middle QPCs in Fig. S11a using G = (R max − R min ) −1 , where R max(min) is the maximum (minimum) value of R at each V bg . The result is shown in Fig. S11c. Note that, since the hBN thicknesses are approximately two times thicker than for Sample 1, the applied V bg and V tg are also larger.
By monitoring R nl = V nl /I vs. B, where V nl is the nonlocal voltage, we obtained the TEF spectra shown in Fig. S11d for V bg = −6 V and at different V tg . This result shows that, even though the background is sensitive to small changes on V tg , the first TEF peak is robust against a V bg = 0.2 V change. Note that R nl shows a peak for B = 0 which is not expected from TEF measurements. The results obtained for V bg > 0 are shown in Fig. S11e and, even though signatures of quantum interference can be seen for B > 0.25 T, no clear signs of TEF can be found.
To confirm that our measurements are in the linear response regime we performed reciprocity checks on the TEF data measured at V bg = −5 V and −6 V. The almost perfect match between both curves confirms that the bias current (100 nA) is small enough so that our measurements are in the linear response regime.
The TEF spectra shown in Figs. S11d and S11f show a peak at B ≈ 0 which would be compatible with the p = 1 TEF peak if the magnet remanence B 0 ≈ −B f . This hypothesis can be ruled out from the reciprocity data. A large B 0 would lead to a horizontal (B) shift of the reciprocal measurements (dashed lines) with respect to the original ones (solid lines). This is a consequence of the reciprocity theorem which states that, in a non-magnetic system and in the linear response regime, R ijkl (B) = R klij (−B), where the first pair of indexes denote the contacts connected to the I source and the second pair the contacts used for the V measurements [18]. If B is shifted by B 0 , then R ijkl will coincide with R klij at B = B 0 instead of 0. Note that both measurements were performed sweeping the magnet from 750 mT to −750 mT so the same B 0 is expected from both measurements. Thus, the coincidence of the direct and reciprocal TEF data at B ≈ 0 confirms that B 0 ≈ 0 and we can conclude that R nl shows an additional peak at B = 0. Most likely, this peak is caused by the detection of the ballistic electron stream reflected at the opposite BLG edge.

S9. NUMERICAL SIMULATIONS
We performed semiclassical calculations of electron focusing in BLG. These calculations require previous knowledge of the angular distribution of the currents and the shape of the Fermi surface. These parameters were obtained using a tight-binding model implemented in Kwant [20]. Both calculations are explained below and the code is available at [5].

A. Tight-binding model
We implement the tight-binding model from [7,21] which includes four hopping parameters: where c n,l and c † n,l are the annihilation and creation operators for electron states at position n and layer l, µ is the chemical potential, ∆ is the layer imbalance, and the sets of hoppings S i , with corresponding strength γ i , are shown in Fig. S12a. Note that γ 2 = 0 in BLG [21] and is not shown in the figure.
From this model, we extract the Fermi surface used in the semiclassical calculations. To reproduce the experimental conditions, we use the displacement field and electron density corresponding to the curves with V bg =+3V in Fig. 1. From the experimental data, we find the corresponding tight-binding parameters µ = 98 meV, and ∆ = 84 meV. We fit the corresponding Fermi surface with the lowest Fourier component that accounts for trigonal warping. Namely, k τ (E, ϕ) = k F,0 + τ δk sin(3ϕ + ϕ c ), where τ = ±1 in valley ±K, ϕ is the polar angle, and k F,0 , δk, and ϕ c are the fitting parameters. We show the computed and fitted data in Fig. S12b.
Finally, because of trigonal warping, it is incorrect to assume that the injected electrons have an uniform angular distribution. Instead, electrons are injected as two valley-polarized jetstreams from the QPC. To obtain the appropriate distribution, we compute the angular distribution of the current density from a QPC with width W i = 7.1 nm using Kwant. The simulated device is depicted in Fig. S12c. We then fit the resulting distribution as where G(θ, θ 0 , Γ) = is a gaussian distribution, θ 0 is the peak position, Γ is the width, and N = π/2 −π/2 (G(θ, θ 0 , Γ) + G(θ, −θ 0 , Γ))dθ the normalization factor. The fitted data is shown in Fig. S12d. Since each gaussian corresponds to one of the valleys, we then use the corresponding sign of theta 0 for each Fermi surface in Eq. S5. We also use a smaller width (Γ/8) to obtain narrow TEF peaks as the ones observed in the measurements.
To obtain the TEF spectra plotted in Figs. 2e and 2f of the main manuscript we assumed that the injector is a point contact with W i ≪ L and the collector has a finite width W c . For every injection angle −π/2 < θ < +π/2, we compute the electron trajectory to determine whether it will reach the detector (that is, for all ⃗ r where y = 0, we determine if L < x < L + W c ). We have assumed that the current dependency on θ follows the distribution shown in Fig. S12d. Each trajectory that hits the collector adds dI/dθ to the collected current. The final result at each B is obtained by summing all the contributions in an equally-spaced distribution of −π/2 < θ < π/2 multiplied by the corresponding dI/dθ.